library(parallel) #并行计算  parApply parLapply parSaplly 
## 利用GenomicFeatures包导入gtf处理
library(GenomicFeatures)

cl <- makeCluster(0.75*detectCores())  #设计启用计算机3/4的核

gene_len <- tibble(gene = rownames(data))

### get gene effective length from gtf ---------
txdb <- makeTxDbFromGFF("~/append-ssd/work/gencode.v43.basic.annotation.gtf",
                        format="gtf")

exons_gene <- exonsBy(txdb, by = "gene") ###提取基因外显子
head(exons_gene)

##计算总外显子长度：用reduce去除掉重叠冗余的部分，,width统计长度，最后计算总长度
exons_gene_lens <- parLapply(cl,exons_gene,function(x){sum(width(reduce(x)))}) 
exons_gene_lens[1:10]

geneid_efflen <- tibble(geneid=names(exons_gene_lens),
                        efflen=as.numeric(exons_gene_lens))

geneid_efflen <- geneid_efflen |>
  mutate(ENSEMBL = str_remove(geneid, '\\..+')) |>
  summarise(efflen = sum(efflen), .by = ENSEMBL)

clusterProfiler::bitr(geneid_efflen$ENSEMBL,
                      fromType = 'ENSEMBL',
                      toType = 'SYMBOL',
                      OrgDb = 'org.Hs.eg.db') |>
  left_join(geneid_efflen) |>
  summarise(efflen = sum(efflen), .by = SYMBOL) |>
  as_tibble() ->
  gtf_symbol_len

write_csv(gtf_symbol_len, '~/work/gencode.v43.gtf-derived-gene-length.csv')

gene_len |>
  left_join(gtf_symbol_len, by = join_by(gene == SYMBOL)) |>
  filter(is.na(efflen))

data[1:5,1:5] |>
  map2(sorted_umi$filter.nUMI[1:5], \(x,y)x/10000*y) |>
  purrr::transpose() |>
  map2(gene_len$askb, \(x,y)x*y)

# use rtracklayer -----------
## slim down gtf for STAR alignment -----------
cellrgr.gtype <- read_csv('00_util_scripts/ref/cellrgr.genetype.csv')

slim_down_gtf <- function(input.gtf, output.gtf){
  rtl.input <- rtracklayer::import(input.gtf)
  
  tb.input <- as_tibble(rtl.input)
  
  allow.gid <- tb.input |>
    filter(gene_type %in% cellrgr.gtype$gene_type,
           transcript_type %in% cellrgr.gtype$gene_type,
           tag != 'readthrough_transcript' | is.na(tag),
           type == 'transcript') |>
    pull(gene_id) |>
    unique()
  
  gene.count <- tb.input$gene_id |>
    unique() |> length()
  
  slim.count <- allow.gid |> length()
  
  message(str_glue('Original gene count: {gene.count}\nFiltered gene count: {slim.count}'))
  
  rtl.input |>
    subset(gene_id %in% allow.gid) |>
    rtracklayer::export(output.gtf)
}

slim_down_gtf('~/append-ssd/work/gencode.vM35.primary_assembly.annotation.gtf',
              '~/append-ssd/work/gencode.m35.pri.reslim.gtf') |> system.time()

# try make a better ENSEMBL:SYMBOL list
rtl_gtf <- rtracklayer::import('~/append-ssd/work/gencode.vM35.primary_assembly.annotation.gtf')

tb.v46 <- as_tibble(rtl_gtf)

tb.v46 |> summarise(n(), .by = transcript_type)
  dplyr::select(gene_id, gene_name, seqnames) |>
  distinct(gene_name, .keep_all = TRUE) |>
  mutate(ensembl = str_remove(gene_id, '\\..+'), .keep = 'unused')

as_tibble(rtl_gtf) |> filter(type == 'gene') |>
  filter(!str_starts(gene_name, '^ENSG'))

write_csv(tb_gtf, 'hg38_gtf-ensembl-symbol.csv')

## T2T version genome!
rtl_gtf.t2t <- rtracklayer::import('~/append-ssd/work/Homo_sapiens-GCA_009914755.4-2022_07-genes.gtf')

vcf.t2t <- as_tibble(rtl_gtf.t2t) 

vcf.t2t |>
  dplyr::select(gene_id, gene_name, seqnames) |>
  distinct(gene_name, .keep_all = TRUE) |>
  mutate(ensembl = str_remove(gene_id, '\\..+'), .keep = 'unused')

as_tibble(rtl_gtf.t2t) |> filter(type == 'gene') |>
  dplyr::select(gene_id, gene_name) |>
  filter(is.na(gene_name))

## cellranger ver gtf
rtl_gtf.10x <- rtracklayer::import('~/append-ssd/work/cellranger_hg38/refdata-gex-GRCh38-2020-A/genes/genes.gtf')

gtf.10x <- as_tibble(rtl_gtf.10x)

gtf.10x |>
  dplyr::select(gene_id, gene_name, seqnames) |>
  distinct(gene_name, .keep_all = TRUE) |>
  mutate(ensembl = str_remove(gene_id, '\\..+'), .keep = 'unused')

gtf.10x |> filter(type == 'gene') |>
  filter(is.na(gene_name))

gtf.10x |>
  summarise(n(), .by = type)
